Part I¶
Variables, if statements, functions, printing output ¶
In [5]:
x<-3
print(x)
In [6]:
if(x<5){
print("The value of x is less than 5")
} else{
print("The value of x is greater than or equal to 5")
}
In [14]:
add5<-function(x){
if(is.numeric(x)){
theanswer<-x+5
return(theanswer)
}
else{
print("Argument x is not a number")
return(0)
}
}
print(paste0("42 plus five is equal to: ", add5(52)))
In [15]:
print(add5("sdafgdsfg"))
In [16]:
running_sum<-0
for (i in 1:10){
running_sum = running_sum+i
}
print(running_sum)
In [17]:
1:10
In [94]:
df<-data.frame(variable1 = c(5,7,9), variable2=c(2,4,6))
df
It's easy to add functions of variables as new variables. For renaming and dropping variables, see the next lesson notebook.¶
In [95]:
df$variable2_squared <- df$variable2^2
df
In [2]:
projectpath<-"C:/Users/Len/Dropbox/Teaching/Data TA/R tutorial/sample project/"
df <- read.csv(paste0(projectpath,"raw data/qcew_wages_industrybyyear_naics.csv"))
In [29]:
df[1:5,c('naics','wage')]
In [6]:
# A "comment" in R can be made with the "#" symbol
# Let's ask R to print how many records are in df:
print(length(df[,1]))
In [10]:
#First five values of wage only:
df$wage[1:5]
In [28]:
print(paste0("Mean industry wage: ", mean(df$wage)))
print(paste0("Max industry wage: ", max(df$wage)))
In [40]:
#What if I want to display these more pretty-like?
print(paste0("Mean industry wage: $", round(mean(df$wage),2)))
print(paste0("Mean industry wage: $", round(mean(df$wage),-3)/1000, "k"))
Suppose we wanted to look up the wage for the coal mining industry in year 2000. The 4-digit "NAICS code" for this industry is 2121. We can do this by slicing df with the logical vector df$naics==2121 & df$year==1990:¶
In [14]:
df[df$naics==2121 & df$year==1990,]
#Be careful, if you used "&&" here it would only apply the AND operation on the first row!
Exercise ¶
What is the average wage for workers in the Logging industry in 2015?¶
Industry codes are listed here¶
Note, this dataset has records for the 4-digit version of industry codes. The longer codes (e.g. 5 digits) are more specific and the data is not dissagregated to that level in this dataset¶
Answer: $41,748.04¶
In [15]:
#Step one: create a dataframe with average manufacturing wage, by year
dfcoal<-df[df$naics==2121,]
dfcoal
We'll use the ggplot plotting library for this, which must be installed. ggplot has a formula based syntax that takes some getting used to.¶
In [16]:
#install.packages("ggplot2")
library(ggplot2)
#ggplot titles by default show up on the left. I like them in the middle:
theme_update(plot.title = element_text(hjust = 0.5))
ggplot(data=dfcoal, aes(x=year, y=wage, group=1)) +
geom_line()+
geom_point()+
labs(title="U.S. Coal Industry Wages Over Time", x="Year", y = "Average nominal wage")
We can also save this plot to a file (e.g. to put it in our paper). Remember the project folder!¶
In [17]:
ggsave(filename = paste0(projectpath,'figures/coal_wagesbyyear.png'), width = 10, height = 5, dpi=500)
In [9]:
projectpath<-"C:/Users/Len/Dropbox/Teaching/Data TA/R tutorial/sample project/"
#split=TRUE allows output to also go to the screen (otherwise you would only see it in the log file)
#append=FALSE tells R to overwrite "mylog.txt" rather than appending to itself, if it already exists
sink(file=paste0(projectpath,'/log files/mylog.txt'),split=TRUE, append=FALSE)
This doesn't work in a Jupyter notebook, so copy this code into R Studio¶
Note: it's good practice to close the connection to this file at the end of your code (this makes more sense when working in RStudio in which your have a ".R" script file). The following code will close any open file connections:¶
In [11]:
sink.reset <- function(){
for(i in seq_len(sink.number())){
sink(NULL)
}
}
Part II¶
Regression analysis ¶
Let's move on to regression analysis. For this, download this data set to your raw data/ folder. Or, you can read it right of the web (and save a local copy):¶
In [15]:
df<- read.csv("http://www.columbia.edu/~ltg2111/teaching/wages_bycounty_yr2000.csv")
write.csv(df, paste0(projectpath,"raw data/wages_bycounty_yr2000.csv"))
This data set contains the average wage and employment count by US county (`fips') in the year 2000:¶
In [9]:
df <- read.csv(paste0(projectpath,"raw data/wages_bycounty_yr2000.csv"))
df[1:5,]
Suppose we want to run a regression of wages on employment in this data. If we just wanted to visualize the regression line, we could let ggplot do it all automatically:¶
In [131]:
#Scatterplot of wage vs employment with regression line
library(ggplot2)
ggplot(data = df, aes(x = employment,y = wage)) +
geom_point(color='cyan') +
geom_smooth(method = "lm", se = FALSE)
Note: the above relationship does not look particularly linear overall. We might want to add employment squared as a regresssor as well, or use so-called non-parametric regression techniques that let the curve take an arbitrary shape.¶
By comparing the linear regression line to the loess local smoother, we can assess the linearity if the relationship. Here we zoom in on counties with average wage between 20-30k and employment below 100k, where the bulk of the data is.¶
In [128]:
ggplot(data = df, aes(x = employment,y = wage)) +
ylim(20000, 30000) + xlim(0,10^5) +
geom_point(color='cyan') +
geom_smooth(method = "lm", se = FALSE) +
geom_smooth(method = "loess", color="green")
Now let's perform the linear regression using the "linear model" or lm command to see the estimated coefficients:¶
In [122]:
lm(wage~employment,data=df)
In [53]:
#It's usually useful to save the estimated model as an object:
model1<-lm(wage~employment,data=df)
summary(model1)
Important! The standard errors that lm computes by default assume homoskedasticity, which is not reasonable in most social science contexts. Based on the above scatterplot, would it be reasonable to assume homoskedastic errors for this regression? Doesn't look like it.¶
To get heteroskedasticity-robust standard errors, you can use the sandwich library: ¶
In [54]:
#install.packages("sandwich")
library(sandwich)
cov <- vcovHC(model1, type = "HC")
model1.robustses <- sqrt(diag(cov))
model1.robustses
In [55]:
model1$coefficients
In [56]:
#If we wanted to supress the constant:
lm(wage~employment+0,data=df)
In [57]:
#Suppose we wanted to inspect the residuals:
hist(model1$residuals, breaks=50,xlim=c(-25000,25000))
In [132]:
temp<- read.csv("https://raw.githubusercontent.com/plotly/datasets/master/fips-unemp-16.csv")
write.csv(temp, paste0(projectpath,"raw data/fips-unemp-16.csv"))
In [133]:
df2<-read.csv(paste0(projectpath,"raw data/fips-unemp-16.csv"))
We can perform the merge using the merge command:¶
In [134]:
merged_data <- merge(df, df2, by="fips", all.x = TRUE, all.y=FALSE)
merged_data[1:5,]
Note that by here can be a list of variables to merge on combinations of, e.g. by=c(var1, var2).¶
The options all.x = TRUE, all.y=FALSE makes sure to keep all of our original county observations, even if we can't find an unemployment number for them. Counties that show up only in the unemployment dataset are ignored.¶
We can see that 9 counties were unmatched, because they have NA's for unemp:¶
In [135]:
summary(merged_data)
merged_data[is.na(merged_data$unemp),]
On the rest of the data, let's do an OLS regression of wage on employment count and unemployment rate (counties with missing unemp will be ignored by the lm command):¶
In [136]:
model2<-lm(wage~employment+unemp,data=merged_data)
summary(model2)
In [137]:
# Again, if we are interested in statistical inference on our coefficients, we should generally use
# standard errors that are valid under heteroskedasticity:
cov <- vcovHC(model2, type = "HC")
model2.robustses <- sqrt(diag(cov))
model2.robustses
In [138]:
head(merged_data)
In [141]:
somevars <- names(merged_data) %in% c("X.x", "X.y")
cleaned_data <- merged_data[!somevars]
head(cleaned_data)
Some other data frame manipulation tasks that may be useful:¶
In [142]:
#Create a new variable that groups employment into 4 quartiles
qnt <- quantile(cleaned_data$employment)
cleaned_data$employment_quartile<-cut(cleaned_data$employment,unique(qnt),include.lowest=TRUE)
#Sort the data frame by employment
cleaned_data <- cleaned_data[order(cleaned_data$employment),]
#Re-order first few columns in data frame.
#In this exmaple, make "employment" and "employment_quartile" first and retain ordering of the rest
somecols <- c("employment", "employment_quartile");
cols <- c(somecols, names(cleaned_data)[-which(names(cleaned_data) %in% somecols)]);
cleaned_data <- cleaned_data[cols]
#rename the "fips" variable to "county"
colnames(cleaned_data)[which(colnames(cleaned_data) == 'fips')[[1]]] <- "county"
head(cleaned_data)
If you have experience with Stata, you'll notice that doing some of these things was more tedious than in Stata.¶
There are some tools that make data manipulation easier in R. One is the dplyr package, another is data.table¶
For example, with data.table, we can quickly extract the average employment across counties within each quartile of employment¶
In [117]:
library(data.table)
dt = data.table(cleaned_data)
dt <- dt[, .(mean(employment)), keyby = .(employment_quartile)]
#the column with mean employment ends up with the non-descriptive name "V1", so let's change it
colnames(dt)[[2]]<-"mean_emp_byqt"
dt
In [118]:
#We may want to merge this information back to the original dataset:
cleaned_data <- merge(cleaned_data, dt, by="employment_quartile")
head(cleaned_data)
Suppose we wanted the average of the variables employment and wage across all counties. We could do this with the lapply funciton, which would allow us to loop through any set of variable names and ask for its mean:¶
In [119]:
nothing<-lapply(c("employment", "wage"), function(thisvar){
print(paste0("Mean of variable ", thisvar, ": ",mean(cleaned_data[,thisvar])))
})
In [23]:
#install.packages("stargazer")
library("stargazer")
stargazer(model1, model1, model2, se=list(NULL, model1.robustses, model2.robustses), column.labels=c("not robust","robust", "robust"), align=TRUE, type="text")
We'll also save a LaTeX version of this table into a file, to use in the next lecture on LaTeX:¶
In [125]:
stargazer(model1, model1, model2, se=list(NULL, model1.robustses, model2.robustses), column.labels=c("not robust","robust", "robust"), align=TRUE, float=F, type="latex", out=paste0(projectpath,"/tables/table1.tex"))
Project!
Suppose we want to add to this regression a "fixed effect" for each state. See if you can figure out how to do this.¶
Some hints:¶
- the first two digits of a FIPS code are constant within a state
- dummies for each value of a variable can be added to the regression by creating a "factor variable", e.g.
df$state.f <- factor(df$state)